Transcriptomic analysis of tuberous root in two sweet potato varieties reveals the important genes and regulatory pathways in tuberous root development

Background Tuberous root formation and development is a complex process in sweet potato, which is regulated by multiple genes and environmental factors. However, the regulatory mechanism of tuberous root development is unclear. Results In this study, the transcriptome of fibrous roots (R0) and tuberous roots in three developmental stages (Rl, R2, R3) were analyzed in two sweet potato varieties, GJS-8 and XGH. A total of 22,914 and 24,446 differentially expressed genes (DEGs) were identified in GJS-8 and XGH respectively, 15,920 differential genes were shared by GJS-8 and XGH. KEGG pathway enrichment analysis showed that the DEGs shared by GJS-8 and XGH were mainly involved in “plant hormone signal transduction” “starch and sucrose metabolism” and “MAPK signal transduction”. Trihelix transcription factor (Tai6.25300) was found to be closely related to tuberous root enlargement by the comprehensive analysis of these DEGs and weighted gene co-expression network analysis (WGCNA). Conclusion A hypothetical model of genetic regulatory network for tuberous root development of sweet potato is proposed, which emphasizes that some specific signal transduction pathways like “plant hormone signal transduction” “Ca2+signal” “MAPK signal transduction” and metabolic processes including “starch and sucrose metabolism” and “cell cycle and cell wall metabolism” are related to tuberous root development in sweet potato. These results provide new insights into the molecular mechanism of tuberous root development in sweet potato. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08670-x.


Introduction
Sweet potato (Ipomoea batatas L) is a dicotyledonous plant of the family Convolvulaceae, growing in tropical, subtropical, and temperate regions, it is the most important rhizome crop after potato and cassava, and one of the most important food crops in the world [1], with an annual global output of more than 100 million tons. China is the largest sweet potato producer in the world, accounting for 80-85% of the global output [1]. Sweet potato is nutritious and contains many ingredients for human health, which has the medicinal values such as anti-cancer, anti-diabetes and anti-inflammatory activity [2,3], and has been selected as one of the test foods for long-term space travel [4]. The tuberous root of sweet Open Access *Correspondence: xiaodong@gxu.edu.cn; lfhe@gxu.edu.cn 1 National Demonstration Center for Experimental Plant Science Education, College of Agriculture, Guangxi University, Nanning 530004, People's Republic of China Full list of author information is available at the end of the article potato is rich in starch and soluble sugar, and its biomass is the highest in all crops. Sweet potato is listed as the key raw material for ethanol production because of its high starch content [5]. How to improve the yield and quality of sweet potato has become a top priority.
Endogenous hormones play an important role in the process of tuberous root expansion. Cytokinin (CTK) and abscisic acid (ABA) are involved in the formation of stored roots [6][7][8][9][10][11], t-zeatin is thought to play an important role in the induction of tuberous roots by activating the primary cambium. ABA regulates the thickening of tuberous roots by activating the cell division of meristem. The content of Auxin (IAA) increased gradually at the initial stage of root expansion in sweet potato tuberous root, and began to decrease after the beginning of secondary growth, while the content of ABA and cytokinin was steadily increased [12,13]. In tuberous root, the content of jasmonic acid (JA) was very high, while the contents in burdock root and fibrous root were less [14].
The growth and expansion of tuberous root in sweet potato are genetically regulated. Previous studies have shown that MADS-box, KNOX genes were highly expressed and related to the expansion of tuberous root in sweet potato [15][16][17]. The overexpression of SRD1 gene promoted the proliferation of cambium cells and xylem cells, and played a role in auxin-mediated initial root thickening [12]. SRF6 was the most abundantly expressed in tuberous root, and its mRNA was located around the primary cambium and meristem of the xylem, promoting the thickening of the tuberous root [18,19]. Besides, an expansin coding gene IbEXP1 was found to play an inhibitory role in the proliferation of cambium cells and xylem cells, which in turn inhibited the initial expansion of tuberous root in sweet potato [19]. The tuberous root development of sweet potato is regulated by multiple genes. However, few genes related to tuberous root development have been identified, and no specific genes regulating tuberous root development of sweet potato have been found, so more researches are needed to reveal the molecular mechanism of tuberous root development of sweet potato.
With the rapid development of sequencing and molecular technology, the study on the molecular mechanism of underlying tuberous root expansion in sweet potato has made great progress. However, the development of tuberous root in sweet potato is a complex biological process, and its mechanism is not clear. Sweet potato is a heterohexaploid plant (2n = 6x = 90) with a genome of 4.4 GB [20]. There are some studies on the development mechanism of sweet potato tuberous root at the transcriptional level. It was found that some specific genes and proteins associated with starch and phytohormone synthesis as well as various transcription factors are involved in storage root formation and development [17,[21][22][23], but there are many genes should be found at transcriptional level. In the meanwhile, previous studies were based on a single variety, however, there are great genetic differences among varieties. It is difficult to explain the general mechanism and variety specificity from transcriptomic analysis using a single variety. In this study, two main sweet potato cultivars with similar developmental processes but having great genetic differences and usually planted in Guangxi Zhuang Autonomous Region of PR China, Xiguahong (XGH, orange flesh sweet potato) and Guijingshu 8 (GJS-8, purple flesh sweet potato), were used as plant materials. RNA sequencing and weighted gene co-expression network analysis (WGCNA) were performed to identify the key candidate genes mediating tuberous root development.

Identification of differentially expressed genes between fibrous root and tuberous root
To explore the molecular mechanism of the formation and development of tuberous roots of sweet potato, 8 cDNA libraries were generated from the fibrous roots(R0) and the tuberous roots at different development stages (R1, R2, R3) in GJS-8 and XGH. Based on Illumina sequencing, a total of 1,514,457,568 original readings were obtained. After removing the connectors, unknown bases and low-quality reads, 1,486,623,198 clean readings were obtained, with an error rate of less than 0.03, Q20 > 97%, Q30 > 93%, which met the quality requirements of database construction. These clean readings were compared to the sweet potato genome using HISAT2 platform, and each library compared the number of reads on the genome to more than 69%. The number of reads aligned to the unique location of the reference genome was more than 63%, and the number of reads aligned to multiple locations of the reference genome was about 3.2-3.8% (Table 1). The sample correlation heat map showed that the R2 value among three biological repetitive samples was greater than 0.8, and that of most of samples was greater than 0.9, indicating that this experiment was highly repeatable and the data were reliable (Fig. 1).

GO and KEGG enrichment analysis of DEGs
To further determine the main biological functions of all DEGs shared by GJS-8 and XGH in the process of tuberous root development, functional annotation was performed by mapping all common DEGs to gene ontology (GO) terms in the GO database. GO enrichment analysis was implemented using a Bonferroni-corrected p ≤ 0.05 as the threshold. Based on this criterion, 33 biological process terms, 3 cellular component terms and 36 molecular function terms were significantly enriched in R1 vs. R0 comparison. Among the DEGs between R1 vs. R0, the "cellular carbohydrate metabolic process" and "single-organism carbohydrate metabolic process" were the major terms of biological process, the "cell wall" and "external encapsulating structure" were the major terms of cellular component, and the "nucleic acid binding transcription factor activity" was the most represented molecular function term (Table S1). A total of 91 biological process terms, 7 cellular component terms and 49 molecular function terms were significantly enriched in R2 vs. R0 comparison. Among the DEGs between R2 vs. R0, the "response to stress" and "single-organism carbohydrate metabolic process" were the major terms of biological process, the "cell periphery", "cell wall" and "external encapsulating structure" were the major terms of cellular component, and the "nucleic acid binding transcription factor activity" was the most represented molecular function term (Table S2). Moreover, 75 biological process terms, 6 cellular component terms, and 39 molecular function terms were significantly enriched in R3 vs. R0 comparisons. Among the DEGs between R3 vs. R0, the "ion transport" and "cell communication" were the major terms of biological process, "cell periphery" "cell wall" and "external encapsulating structure" were major terms of cellular component, and "nucleic acid binding transcription factor activity" was the most represented molecular function term (Table S3).

Comprehensive analysis of differential expression of signal transduction pathway genes
The KEGG enrichment analysis of the DEGs shared by GJS-8 and XGH during tuberous root expansion showed that they were significantly enriched in many signal transduction pathways. Furthermore, these DEGs were annotated using NR, GO, and KEGG annotations, Pearson correlation between samples

Fig. 1 Correlation analysis between sample replicates
and a large number of DEGs were involved in signal transduction, cell wall, cell division, starch and sucrose metabolism pathways, indicating that signal transduction pathways played an important role in the process of sweet potato tuberous root expansion. Therefore, we analyzed the related genes of these pathways.

Hormone signal
In this study, a total of 58 genes related to biosynthesis, metabolism and signal transduction of various hormones were identified (Table S4). The auxin signal transduction pathway was the most active, followed by ethylene signal transduction pathway. The genes related to hormone signal transduction in two varieties at the same developmental stage were further analyzed. In the auxin pathway, 3 The number of differentially expressed genes between tuberous root and fiber root at different stages of GJS-8 and XGH. a Statistics on the number of differential genes in different situations. b The number of differential genes between GJS-8 and XGH calmodulin-binding protein (CaM/CaM-binding), and 26 Calreticulin (CBL) were identified from the common DEGs of two varieties (Table S6). It is worth noting that most of genes were down-regulated in whole expansion stage. A total of 22 phospholipid signal-related genes were identified from the common DEGs of two varieties (Table  S7). Among them, 6 genes were significantly up-regulated in R1, R2 and R3 stages, 4 genes were significantly downregulated in R1, R2 and R3 stages.

Cell wall and cell cycle
We identified 95 genes related to cell wall and cell cycle from the DEGs shared by GJS-8 and XGH (Table S9)  kinases (CDKs) and 11 cyclin-dependent kinase inhibitors (CDKIs). Among these genes, most of XTH and CDC genes were down-regulated, and most of the genes related to FtsZ, CDC5 and CDKIs were up-regulated.

Starch and sucrose metabolism
Seventy genes related to starch and sucrose metabolism were identified from the DEGs shared by GJS-8 and XGH (Table S10)

Transcription factor
In this study, 296 TF genes were identified as DEGs shared by GJS-8 and XGH. Among them, 126 TFs were up-regulated, and 170 TFs were down-regulated during the tuberous root development (Table S11). WRKYs, HBs, MYBs were the major represented TF families (Fig. 4). Twenty-nine transcription factors in these families were significantly up-regulated, and their expression levels increased successively in the R1, R2 and R3 stages of tuberous root development in two cultivars, it mainly included the family of HB, C2H2, MYB transcription factors (Fig. 5).

Weighted gene co-expression network analysis
To further understand the relationship between gene expression and tuberous root development, the weighted gene co-expression network analysis (WGCNA) was performed. In this study, β (soft-power threshold) = 9 was set to guarantee high scale independence and low mean connectivity (near 0) (Fig. 6A). The dissimilarity of the modules was set as 0.75, and a total of 14 modules were generated (Fig. 6B). The module trait relationship was shown in Fig. 6C. Green modules are highly related to tuberous root development (r > 0.80, p < 0.005). GO enrichment analysis was further carried out on the genes of green module (Table S12). The result showed that the biological processes were the most enriched in this module related to energy metabolism and transport. In addition, it was also significantly enriched in mRNA processing, hormone response, endogenous stimulus response and stress response. KEGG enrichment analysis showed that the green module was significantly enriched in transcription factors, plant circadian rhythm (sot04712), MAPK plant signal pathway (sot04016), and plant hormone signal transduction (sot04075) ( Table  S13). The gene connectivity in the modules represents the regulatory relationship between the gene and other genes. The higher the connectivity, the greater the regulatory role of the gene in the modules, the more likely it was a hub gene. The gene with the highest connectivity in the green module was selected as the core gene of the module. This gene encoded a trihelix transcription factor (Tai6.25300). The homology of this gene in Arabidopsis is AT1G13450.1 (trihelix transcription factor: GT-1). A total of 1272 genes interacted with trihelix, including genes related to light signaling, calcium signaling, and plant hormone signaling, implying the processes the genes involved were potentially co-regulated. The interaction network of core genes was visualized by Cytoscape software. Because there were many genes interacting with hub genes, only partial genes were shown here (Fig. 7).

Genes with significant differences in tuberous root development between two varieties
Taking | log2 (FoldChange) | > 1 and padj < 0.05 as the standard, we identified 18,028 differentially expressed genes (DEGs) for the GJS_8 vs. XGH (R1, R2 and R3), of which 12,792 were in R1 stage, 9979 in R2 stage and 8828 DEGs in R3. KEGG enrichment analysis showed that the up-regulated genes were significantly enriched to phenylpropanoid biosynthesis (sot00940), flavonoid  (Table S14). The difference of these anthocyanin related genes was the greatest in the R1 stage of the two varieties, and the difference was more than 10 times.

Verification of gene expression patterns by qRT-PCR
In order to verify the accuracy of RNA-Seq results, we randomly selected 6 genes (Tai6.25300, Tai6.22648, Tai6.3107, Tai6.42353, Tai6.46822, and Tai6.24971) for qRT-PCR analysis. The results showed that the expression pattern of these 6 differential genes was similar to that of RNA-Seq (Fig. 8). The results indicated that the RNA-Seq was reliable.

Discussion
The formation and development of the tuberous root of sweet potato is a complex process, which mainly involves the formation of vascular cambium and secondary cambium. After the formation of round vascular cambium, the tuberous root begins to thicken, then the cells continue to proliferate and expand to form a secondary cambium, which is accompanied by the continuous accumulation of starch and other substances, resulting in the continuous enlargement of the tuberous root. Previous studies showed that the meristems are always active during tuberous root bulking, the transcriptome data obtained in this study reveal that the regulators of meristem development, such as LBD4 (LOB domaincontaining protein 4, Tai6.18322, and Tai6.27010), WOX4 (WUSCHEL HOMEOBOX RELATED 14, Tai6.17770, and Tai6.44989) were significantly upregulated at tuberous root development, which is consistent with the results of previous studies [23]. Moreover, the genes are involved in cell division, including cell division protein FtsZ (FtsZ), cell division cycle 5 (CDC5), cell division control protein (CDC), and cyclin-dependent kinase (CDK), their expression levels were significantly enhanced in the tuberous root expansion stage (Table  S9). The Genes involved in cell extension and expansion, including extension, XET, and expansin, also were significantly enhanced in the tuberous root expansion stage (Table S9). These results indicate that the formation and development of tuberous roots are inseparable from the active meristems and cell division.
A series of studies have shown that the initiation and induction of root/tuber is affected by the environment. For potatoes, photoperiod is essential for tuber formation [24]. Moreover, light is also important for the expansion of Rehmannia glutinosa tuberous root [25]. Photoperiod response protein, lateral organ boundaries protein (LOB), and GATA transcription factor are important members of photoperiod regulation. In this study, the expression of LOB (Tai6.27900) and GATA (Tai6.27468) were significantly enhanced during the tuberous root expansion stage. Furthermore, genes related to light signal transduction including phototropin, CONSTANS, and COP-interactive proteins were also significantly enhanced during the tuberous root expansion stage (Table S8). However, their peaks and expression patterns were obviously different, suggesting that light regulation is very critical to tuberous root formation and continuous development.
Moreover, genes detected in the roots may also be transcribed in the leaves and then transported to the root. For example, after being transcribed in leaves, potato stBEL5 mRNA was transported through the phloem to the stolon tip for translation into protein, thereby promoting the formation of storage organs [26]. In this study, 14 BELs genes were consistently up-regulated during the tuberous root expansion stage (Table S8), which suggest that these genes may be functionally similar to the stBEL5. Although the storage organs of potato and sweet potato are different, they may have similar regulatory systems. Therefore, they may be involved in light signal-regulated tuberous root development via similar mechanisms.

The relationship between hormones and tuberous root swelling
Hormones are important signals in plant root development [27,28]. In this study, the plant hormone signal transduction pathway was one of the most enriched KEGG pathways in tuberous root expansion stage. Auxin plays an important role in cambium cell proliferation and cell expansion [12], also maintains the meristem state of cambium cells and increase the number of xylem elements [29]. In the studies of radish, Rehmannia glutinosa and Callerya speciosa, the expressions of auxin-related genes were significantly up-regulated during tuberous root expansion stage [25,30,31]. In this study, 7 auxinrelated genes (AUX / IAA, ARF, SAUR, and CH3) were up-regulated in tuberous root expansion stage, implying that they may relate to cell expansion in the secondary growth of cambium.
The results showed that cytokinin was involved in the proliferation and development of cambium cells, and the expression reached the highest level in the rapid growth stage of tuberous root, which was related to the development and formation of tuberous root / tuber [29,[32][33][34]. In this study, the expression of cytokinin related gene (Tai6.10485) was significantly up-regulated during tuberous root expansion, suggesting that cytokinin may promote root expansion by participating in the development of cambium.
Ethylene is a key regulator of rhizome induction and development [35], which promotes tuber formation by inhibiting GA biosynthesis [36]. Moreover, it has been shown that GA, auxin,and ethylene affect cell growth in the root by opposing the action of DELLA proteins. In this study, the expressions of ethylene-related genes were significantly up-regulated during tuberous root expansion (Table S4). Overall, these results suggest that these hormone signals related genes play vital roles during the tuberous root expansion stage.

Multiple signal pathways are activated to regulate tuberous root development
Cellular processes involved in a series of signaling pathways are usually triggered by specific stimuli and hormones. Phospholipid signal plays an important role in root growth, cell division, and hormone regulation [37,38]. It was reported that the expression levels of phospholipid signal-related genes/proteins were increased in the early stage of tuberous root expansion in Rehmannia glutinosa. In addition, the phospholipid-calcium signal system regulated potato tuber formation [25,39]. In this study, 6 phospholipid signal-related genes were up-regulated in the stage of tuberous root expansion in GJS-8 and XGH, and the expression profiles in two varieties were quite similar, indicating that phospholipid signal was involved in the initiation and of tuberous root expansion.
Calcium is one of the main nutrients and is involved in almost the whole process of plant growth, including the controls of cell division, differentiation, and stress response as the second messenger [40,41]. Studies revealed that CDPK played a role in the signal pathway of root initiation in potato and cassava, and exogenous calcium levels could affect the quantity and weight of potato tuber [42][43][44]. In addition, Ca 2+ concentration and calcium signal-related genes (CBP, CBL, CaM, and CDPK) were significantly up-regulated during tuberous root formation in Rehmannia Glutinosa [25]. In this study, there was an increase in the stage of tuberous root expansion in the expression level of calcium signalingrelated genes, including 9 CDPKs, 8 CBLs, and 1 CaM (Table S6), which suggests that calcium signal is involved in the formation and expansion of tuberous root in sweet potato. In addition, some genes related to the MAPK signaling pathway were up-regulated during tuberous root expansion development (Table S5), suggesting that the MAPK signal participats in the initiation and expansion of tuberous root formation. It has been shown that the MAPK signal plays an important role in cell cycle regulation, hormone, and stress response [45].

Transcription factor regulation and weighted gene co-expression network analysis
Transcription factors play an important role in the regulation of plant growth and development and secondary metabolism. Many transcription factors have been identified to play key roles in organ development, including MADS, bHLH, MYB, NAC, GRAS et al. In this study, we identified 29 transcription factors that were significantly up-regulated during the tuberous root expansion stage in two varieties. Their expression levels increased successively (Fig. 5). Among these TFs, MYBs and HBs were the main transcription factors with large up-regulation multiples. One trihelix transcription factor gene (Tai6.25300) was identified as a tuberous root expansion-related gene through WGCNA analysis, its homologous gene in Arabidopsis was AT1G13450.1(Trihelix, GT-1), which was considered to be a molecular switch responded to light signals through Ca 2+ -dependent phosphorylation/ dephosphorylation [46]. The Trihelix factor is a plantspecific triple helix DNA binding transcription factor. Many studies have proved that the trihelix transcription factor was involved in plant light response [47,48]. In this study, the expression of light signal related-genes was coordinated with Tai6.25300, and significantly up-regulated during tuberous root development. Moreover, qRT-PCR confirmed that the expression of Tai6.25300 was up-regulated and increased successively during tuberous root development in two varieties, suggesting that Tai6.23500 was closely related to tuberous root development. We infer that Tai6.25300 participates in tuberous root expansion by positively regulating light signal related genes.
MYBs were involved in cell cycle regulation, plant morphogenesis, cell wall synthesis, secondary metabolism, xylem/phloem differentiation, root radial pattern formation, and so on [49,50]. Furthermore, previous studies have found that the transcriptional level of MYB was significantly up-regulated during rhizome development [30,51], and MYBs were highly expressed at the rapid thickening stages of Callerya speciosa [36]. In this study, 32 MYB transcription factors were significantly differentially expressed as tuberous root development, of which 8 were significantly up-regulated. Homeodomain (Homebox, HB) transcription factors are very important regulatory proteins in plants, which are mainly divided into 14 categories, including KNOX, BEL, and HD-ZIP, etc. Arabidopsis HB transcription factors were involved in cell division, differentiation, replication, growth, and regulation of the early development of vascular tissue [52,53], In addition, the members of HB family were also involved in the regulation of cambium cell differentiation to phloem and lignin biosynthesis [54,55]. RNA-Seq data revealed that 3 homeobox genes were notably upregulated during the formation and thickening of storage roots [22]. In this study, 36 HB transcription factors were significantly differentially expressed in tuberous root development, of which 26 were significantly up-regulated.
To sum up, these results suggest that transcription factors may drive root/stem growth through cell cycle regulation, cell division, and secondary wall strength. The TFs revealed in this study may be the important candidate genes for breeding sweet potato with high production in the future.

Starch and sucrose metabolism regulation
Sucrose and starch accumulation occurs during the bulking of storage roots, they are considered to be one of the most important carbohydrates, and play an important role in the formation of storage organs. Sucrose invertase and sucrose synthase were involved in the introduction and accumulation of sucrose in storage roots [56]. In addition, sucrose synthase was related to the tuber / tuberous root growth of potato and radish and was a key enzyme in the early development of radish storage root [57][58][59][60]. In this study, 5 SuSy genes were significantly upregulated during tuberous root development in GJS-8 and XGH, while 2 INV genes were significantly downregulated (Table S10), Invertase was active in fibrous roots of sweet potato but rapidly decreased to an undetectable level during storage root development [61]. Furthermore, Jackson showed that high content of sucrose was required as a necessary condition during the formation of storage organs [62]. In the present study, SPS (Tai6.24187), the major source of sucrose synthesis activity [63], was up-regulated during tuberous roots expansion. This result was consistent with previous studies in radish that found up-regulation of SPS playing a major role in the thickening stage of radish taproot [64].
The accumulation of starch occurs at the same time as the expansion of storage organs. It has shown that the expansions of potato and lotus root tubers were highly coordinated with the accumulation of starch [65,66]. The expansion of cassava root was synchronized with the accumulation of starch [67], and granule-bound starch synthase (GBSS) has been shown to affect starch synthesis in storage organs [68]. In this study, 22 starch-related genes (6 GBSSs, 4 SSSs, 8 SBEs, and 4 isoamylases) were significantly up-regulated during root tuber expansion (Table S10), which was similar to previous studies. SBE, GBSS, and SS-related genes were significantly up-regulated during root expansion of Panax notoginseng [69]. These starch and sucrose metabolism genes play important roles in tuberous root expansion.

Genes with significant differences in tuberous root development between two varieties
GJS_8 and XGH are two varieties with different anthocyanin content. GJS_8 has higher anthocyanin content than XGH. Anthocyanins are water-soluble pigments and an important class of flavonoids. We found that there was a large number of genes with significant differences in tuberous root development between two varieties. KEGG enrichment analysis showed that the DEGs were significantly enriched to phenylpropanoid biosynthesis (sot00940), flavonoid biosynthesis (sot00941), and starch and sucrose metabolism pathway (sot00500). It was also found that phenylpropanoid biosynthesis and flavonoid biosynthesis was significantly enriched in the process of anthocyanin biosynthesis [70]. In addition, we identified a large number of MYB, bHLH, WD40 transcription factors, and anthocyanin biosynthesis genes from these differential genes, including 6 MYBs,17 bHLHs, 3 C4Hs, 5 4CLs,6 CHSs, 2CHILs, and 2 LDOX/ANSs, which were significantly differentially expressed between GJS_8 and XGH and also significant differentially expressed between tuberous root and fiber root, especially in GJS_8 tuberous root. A large number of studies have shown that MYB, bHLH, and WD40 transcription factors were the regulators of flavonoid biosynthesis, and the results also showed that IbMYB1 controls the biosynthesis of anthocyanins in sweet potato [71]. It was found that 10 anthocyanin biosynthesis genes were significantly up-regulated during Aronia melanocarpa fruit development [72]. Hence, it shows that anthocyanin biosynthesis related-genes may be involved in the tuberous root development in sweet potato, and their regulatory mechanism should be studied in the next step.

Regulatory networks associated with tuberous root development
Tuberous root development is a complex regulatory process, which is affected by many factors. In this study, through transcriptome analysis, combined with previous research results, a hypothetical model of sweet potato tuberous root development regulatory network is proposed (Fig. 9). The cells in the vascular cambium divide and expand continuously to produce secondary xylem and secondary phloem, resulting in the expansion of tuberous root. Cell proliferation is regulated through several signal transduction pathways (light, Phospholipid, calcium, MAPK, hormone, and transcription signaling) and metabolism possesses (cell wall, sucrose, and starch metabolism). Several genes including photoperiod (LOB, GATA, Phototropin, COL, and COP), calcium signal (CDPK, CBL, and CaM), MAPK signal, auxin-related genes (Aux/IAA, CH3, ARF, and SAUR), HB transcription factors (BELL, KNOX, and HD-ZIP), are highly expressed to promote cell differentiation, division, expansion and sucrose and starch accumulation at the secondary structure. In addition, FtsZ, CDC, CDK, XTH, expansin, and extension, are involved in cell division extension and expansion. Finally, SuSy, SPS, SSS, GBSS, and SBE are involved in the hydrolysis of sucrose and the synthesis of starch. Further functional identification studies were needed to confirm the functions of these potential genes.

Conclusion
Integrated transcriptomic and WGCNA analyses were performed in the study, there were 15,920 differential genes shared by XGH and GJS-8. GO and KEGG pathway enrichment analysis revealed that these DEGs were mainly involved in plant hormone signal transduction, starch and sucrose metabolism, MAPK signal transduction, light signal, phospholipid signal, calcium signal, transcription factor, cell wall, and cell cycle. Furthermore, WGCNA and qRT-PCR analysis suggested that Tai6.25300 played an important role in tuberous root development in sweet potato. A hypothetical model of a genetic regulatory network associated with tuberous roots in sweet potato is put forward. The tuberous root development of sweet potato is mainly attributed to cell differentiation, division, and expansion, which are regulated and promoted by certain specific signal transduction pathways and metabolism processes. These findings can not only provide novel insights into the molecular regulation mechanism of tuberous root expansion, but also support theoretical basis for genetic improvement of sweet potato.

Materials
Two sweet potato varieties, GJS-8 and XGH were used in this study. They were planted in the experimental farm of Hepu Institute of Agricultural Science in Beihai, Guangxi. At 90 days after planting, Sample collection refers to Ku et al's method [14], Fibrous roots (R0:RGJ8_0, RXGH_0; 1 mm diameter) and developing tuberous roots [(R1:RGJ8_1, RXGH_1; 1 cm diameter, less than 2 g), (R2:RGJ8_2, RXGH_2; 3 cm diameter, 5-10 g), (R3:RGJ8_3, RXGH_3; 5 cm diameter, approx 50 g)] were collected,respectively (Fig. 10). Three plants were selected randomly from every repetition each time. At least five roots were mixed as a biological biological repetition. For the big tuberous root samples, five fresh tuberous roots from a repetition were washed with distilled water, cut down into slices, and mixed as a biological repetition. Three biological replicates were performed. The samples were stored at − 80 °C for extracting total RNA.

RNA extraction, cDNA library construction, and RNA-Seq
A conventional trizol method was used to extract RNA from the samples. The concentration and purity of total RNA were determined by a NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (AgilentTechnologies, CA, USA). Sequencing libraries were generated using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, USA).

RNA sequencing and data analysis
3 μg total RNA from each sample was used as the input material, fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/ polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 250 ~ 300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High -Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. Clean reads were obtained by removing reads containing an adapter, reads containing ploy-N and low-quality reads from the raw data. The clean reads were then aligned with the sweet potato genome (http:// public-genom es-ngs. molgen. mpg. de/ cgi-bin/ hgGat eway? hgsid= 9052& clade= plant & org= Ipomo ea+ batat as& db= ipoBa t4 ) [23]. Feature Counts v1.5.0-p3 was used to count the read numbers mapped to each gene, and the FPKM of each gene was then calculated based on the length of the gene and the read count mapped to the gene [23]. Genes with an adjusted P-value < 0.05 and | log2 (FoldChange) | > 1 obtained by DESeq2 were considered DEGs.

Functional annotation
Gene Ontology (GO) enrichment analysis of the DEGs was implemented using the cluster Profiler R package, and the gene length bias was corrected during this process [73]. KOBAS software was used to test the statistical enrichment of the DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [74]. To obtain more information about the DEGs, the DEGs were annotated using seven databases: NR (NCBI nonredundant protein), NT (NCBI Nucleotide Sequences), Gene Ontology (GO), KO (KO, KEGG Orthology), KOG (Eukaryotic Or Thologous Groups), Pfam (Protein Family Database) and Swiss-Prot (a manually annotated and reviewed protein sequence database). All the DEGs were subjected to hierarchical clustering analysis using the average linkage method [75].

Weighted gene co-expression network analysis
The DEGs detected with DESeq2 were combined and the TPM values for the 24 samples were determined. Each TPM value was increased by 0.01 and further transformed by a log10 calculation. The converted data were analyzed with the R package WGCNA (version 1.66), with a power value of 9 [76,77].

Validation of the DEGs data using qRT-PCR
Total RNAs were extracted from the tuberous samples (fibrous root, tuberous roots less than 2 g, tuberous roots 5-10 g, tuberous roots greater than 50 g) with Trizol ® Reagent (Magen, China). and then reverse transcribed into cDNA with HiScript III SuperMix for qPCR(+gDNA wiper) (Vazyme, China). qRT-PCR was carried out using SYBR Premix Ex TaqII Kit (TaKaRa, Dalian, China) on a Bio-Rad iQ5 Real-time PCR System (Bio-Rad Laboratories, CA, USA), Ten μl reaction solution contained 5 μl SYBR Green I Master, 1 μl specific Primer, 1 μl cDNA samples, 3 μl RNase-Free H2O. One-third dilution of the cDNA sample was used, and the reaction conditions were: 30s at 95 °C followed by 40 cycles of 30s at 95 °C, and 30s at 60 °C. Each sample had three biological replicates with three technical replicates for each biological replicate. The relative expression level was calculated by the equation ratio 2 -ΔΔCt . The primers of selected genes were designed using primer 5 software (Table S15), and UBI gene was used as the internal control.